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Abstract 

Within the framework of time-dependent density functional theory combined with the Korringa- 
Kohn-Rostoker Green function formalism, we present a real space methodology to investigate dynam- 
ical magnetic excitations from first-principles. We set forth a scheme which enables one to deduce the 
correct effective Coulomb potential needed to preserve the spin-invariance signature in the dynamical 
susceptibilities, i.e. the Goldstone mode. We use our approach to explore the spin dynamics of 3d 
adatoms and different dimers deposited on a Cu(OOl) with emphasis on their decay to particle-hole 
pairs. 
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I. INTRODUCTION 



The magnetic functionalization of nanostructures made of few atoms requires the un- 
derstanding of spin-excitations at the nanoscale and subnanoscale level. Recently, state of 
the art experiments based on scanning tunneling microscopy (STM) were utilized to excite 
and control the magnetic states of single adatoms sitting on semi-insulating [l| or metal- 
lic 2], Q] surfaces. The spin dynamics of moment bearing 3d metal atoms have been probed 
in those experiments but often the theoretical picture used for the interpretation is based on 
a model Hamiltonian describing an atomic like localized moment with integer or half-integer 
spin. Such a model is useful only for systems where the substrate interacts weakly with the 
adsorbate it fails qualitatively to describe cases with strong coupling to the substrate 
electrons where hybridization leads to moments far from integer and half integer values, and 
d levels with widths that can range from a few hundred millivolts to perhaps an electron volt. 
This paper presents a scheme wherein one may address the commonly encountered strongly 
coupled systems, with density functional theory as the basis. In contrast to empirical tight 
binding schemes used earlier I4- 7| the method set forth in this paper incorporates a proper 
ab-initio based description of the one electron physics from upon which our description of 
spin dynamics is erected. Also the scheme set forth in this paper may be implemented with 
modest computational labor. 

Several approaches have been proposed to describe inelastic STM experiments involving 
the above mentioned local moment picture ls-12] but none are based on taking full account 
of the electronic structure of the adsorbates as well as the substrates including the effects 
of hybridization. The latter requires, among other ingredients 13], the evaluation of the 



transverse magnetic response function x or the* so-called transverse dynamical magnetic 
susceptibility that relates, in linear response theory, the amplitude of the transverse spin 
motion m XiY produced by a transverse external magnetic field B ext of frequency u. There 



are three major roads followed to compute X- 0) empirical tic 
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(ii) time-dependent density functional theory (TD-DFT) 



ht-binding theory (ETB) 



EL 



14|-|21I], and (iii) many b od\ 



perturbation theory (MBPT) using the Random Phase Approximation (RPA) and DFT [21 



231 ]. The calculation of x requires one to solve a Dyson equation whose solution may be 
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written schematically in the form: 

X = Xo(l - ^Xo)" 1 (1) 



As noted in Ref. [211 ], %o is described by a di 



used to calculate it. Within TD-DFT 
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ferent nomenclature depending on the method 



151 ], Xo is known as the Kohn-Sham susceptibility 



and U is the exchange and correlation kernel that if ideally known completely would render 
Eq. [T]the exact solution. U is obviously approximated in practice, for example, by the adia- 
batic local spin density approximation. It turns out that evaluating Eq.[T]is computationally 
very challenging, especially within the TD-DFT or the MBPT. This explains the very few 
calculations found in the literature, almost all of which address bulk systems. This makes 
it even more challenging to simulate inelastic STM experiments that examine adatoms de- 
posited on surfaces. Recently, we developed a method [21] that handles the calculation of the 
transverse dynamical magnetic susceptibility in a scheme that resembles ETB but is based 
on TD-DFT. Thus the method incorporates full self-consistent first-principles calculations of 
the underlying electronic structure. Two interesting results were obtained: (i) a justification 



of the Lowde and Windsor scheme 24] emerged from the analysis and (ii) values of U deter- 
mined from first-principles for different systems are in accordance with the empirical values 
extracted from photoemission data by Himpsel. 
In our previous paper 



211 ], we addressed a central question related to the practical deter- 
mination of Xo and U, within the framework of density functional based schemes. It is known, 
but often not discussed explicitly, that the Goldstone theorem is not satisfied, in practice, 
when solving Eq. [1] within TD DFT schemes. We remark that within the framework of the 
empirical tight binding scheme, the Goldstone theorem is satisfied exactly, as demonstrated 
earlier Q]. The Goldstone theorem, when satisfied, insures that the zero wave vector spin 
waves have precisely zero frequency (when spin orbit coupling is set aside). The reason the 
Goldstone theorem is not satisfied within density functional based schemes is that the nu- 
merical methods used to extract U and Xo are not compatible with the Ward identity. To 
compensate for this problem, Sasioglu et al. [23] correct U by 45% in their study of bulk 
Ni while Buczek et al. find a finite frequency for the Goldstone mode [19]. To cure such 
inconsistencies, an ad-hoc shift by hand of the value of U is used commonly. Our aim is to 



demonstrate that such corrections could be dangerous, for instance, when the system under 
investigation contains more than two atoms in the unit cell. In Ref. [21|] , we set forth and 



utilized a sum rule that allows one to generate a U that is fully compatible with the Goldstone 
mode. 

The discussion of the sum rule in ref. [21|] was brief, though its application was illustrated. 



2 II ] including the sum rule needed 



In this paper, we provide a detailed derivation of our scheme 
to determine U. Our method is based on the Korringa-Kohn-Rostoker single particle Green 
function (KKR-GF) 



26| which contains an ab-initio description of the electronic structure. 
We remark that in earlier work, the empirical tight binding method has been used suc- 
cessfully to describe spin waves in films on substrates_ 5, Q] along with the spin dynamics of 
adatoms as probed in the recent STM experiments In this approach, it is necessary to 
make contact with electronic structure calculations for the purpose of extracting the tight 
binding parameters required to describe the one electron properties of the system of interest. 
Often appropriate electronic structure calculations are unavailable, or if they are it can be a 
challenge to extract appropriate parameters in an unambiguous manner for complex systems 
such as ultrathin films adsorbed on substrates. The approach we develop here eliminates 
this issue completely, while at the same time it provides a computationally straightforward 
scheme for generating the dynamic transverse susceptibility. 



II. STRUCTURE OF THE THEORY; THE SUM RULE AND THE EFFECTIVE U 

It is, of course, possible in principle to calculate the Kohn-Sham non interacting suscepti- 
bility Xo- l n this section, we show that once Xo assumed known, we can derive a prescription 
for generating the effective Coulomb interaction U which enters Eq. [T]that is fully compatible 
with the Goldstone theorem. In effect, U is a functional of Xo- With U determined in the 
manner we describe, there is no reason for ad-hoc adjustment of this central parameter. We 
also describe a scheme which allows one to generate a physically sensible approximation to 
Xo that is straightforward and simple to implement. We then use this scheme to gener- 
ate a series of explicit predictions regarding the nature of spin excitations of adatoms and 
adatom-dimers. 
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To begin, we assume we have in hand a magnetic system with an initial charge density 
n Q (r). Its ground state magnetization (m z (F)) pointing along, say, the z-direction experi- 
ences a modification induced by a small time-dependent external transverse magnetic field 
S ex t(r; £). The result is an induced transverse magnetization m xy (r; t) localized in the (xy) 
plane perpendicular to the direction z. To describe the transverse magnetization, we begin 
by calculating the frequency dependent Kohn-Sham transverse susceptibility or Xo which may 
be expressed in the form 

xUr,f;u) = ~J dzfiz^Gfjif.f-.z + vj^mGjiif.rz) 

+ \mG\ j {r,7;z)GJ i x {7,r,z-u)) (2) 

where f (z) is the Fermi distribution function, G and G~ represent the retarded and advanced 
one particle Green functions connecting atomic sites / and j and ImG = — ^(G — G"). 

A comment on the notation is in order. In general, the point r* is in unit cell /', and 7 is in 
unit cell j. These vectors are measured from the center of their respective unit cells. Thus, 
if we wish to describe these two points with respect to a master origin O, we will describe 
the notation f+Rj and 7 + Rj, respectively where R jd are vectors from O to the center 
points of cell /', j. With this convention in mind, the single particle Green function, often 
described as G(r + R h 7 + Rj,z), will here be described as G (J (r, f\z), a notation that is 
very convenient when the KKR scheme we employ is utilized. 

To derive our criterion for choosing an effective U, our interest is in the static form of 
the Kohn-Sham susceptibility. At u = 0, the expression in Eq. E] reduces to the usual form 
of the static magnetic susceptibility: 

7&{T.r;Q) = dzf(z)(Gf J (r,r t ;z)Gj l (r,r,z) 

- Gr j \r,f;z)G^{f,r,z)) (3) 

Our first step it to multiply both sides of Eq. |3]by B J eff (F;u = 0) and then we integrate 
over f within the atomic site j and sum up over all sites j: 

£ / ^(^0)^(^0) = ■£-[ dzf(z)J2 I df (4) 

j j 

(Gl.(r,?;z)Bi ff (7;0)G} i (r',r,z) 
- G^(r, 7; z)B{ ff {7; 0)G~\7, r; z)) (5) 
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B eff is given by the difference between the potentials of each spin channel ( V± — Vf). 

We next use an identity derived in the Appendix that relates the Green function for a given 
spin channel, say f, to the Green function of the opposite spin channel through the potential 
difference B eff : 

Gl(r, n z) = Gj(r, r; z) + W drG^K 7; z)B{ ff {f- 0)Gj (r<, r; z) (6) 

j J 

or 

Gl(r, n z) - Gf,(r, r; z) = £ / ^Gj(r, f : z^f; 0)Gt(r<, r; z) (7) 

j J 

Similar relations but written differently have been already used for example in Refs. 
Thus Eq. [5] becomes: 

(6j(r.?;z)-6i(r.^z) 

- G^^zJ + G^Cf.^z)) (8) 
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which is the same as 



d?x!j(f,?; 0)Bi ff (f;0) 



(lmGt(r,r;z)-lmG, J ;.(r,r;z)) 



(9) 



One can recognize that the right-hand side of the previous equation is simply m' z (F; 0). Thus, 
we obtain the final form of an important sum rule: 



£ y G^'(r,r ;W = 0)6^(^ = 0) = m> z {r,u = V) 



(10) 



We remark that within the empirical tight-binding scheme, a statement equivalent to 
EqnU]is found in Ref. f|. 

The Kohn-Sham susceptibility x'o(rt, can be expanded in terms of real spherical 
harmonics, Y and when this is done it can be expressed as a sum over angular momenta L, L\, 
L 2 and L 3 as Zlul 2 l 3 X^ LLl ' JLzL!l {rt. r' t')Y L {f)Y Ll {?')Y L2 {r)Y L3 {?) . This follows since X o is 
a convolution of single particle Green functions (see Eq.ED- Consequently, within the atomic 



sphere approximation (ASA) and assuming a spherical magnetic field m'(?t) = m' xy (rt), 
m' z (rt) = m' z (rt) and B J ext (Pt) = B J ext (r't), Eq. [TO] reads: 

Hl dr ' E Vl(OVl3(Ox^ i;jL2L3 (r,r';0)ei ff (r';0) x 
j J LL1L2L3 

J dfY Ll (f)Y L2 (f) = m^r.O) (11) 

If one integrates both sides of the previous equation over dr and uses J d?Y L (?)Yy{?) = 5i_u 
one finds: 

E / dr'^X l o LLljLlL (r,r':0)Bi ff (r':0) = 47r<(r;0) (12) 

If we define 

= flU(^o) (13) 

4W z (r';0) 

that is the usual form for the effective U that enters Eq. [I] as generated from the Adiabatic 
Local Spin Density Approximation given in the upcoming section, then Eq.[T2]can be rewritten 
as: 

j dr'^Xo LLi;jXlL (^0)m^ z (^0)L/ J (r') = <(r; 0) (14) 



j " LL 



or as 



[ dr'r u (r,r')U J (r') = m' z {r; 0) (15) 

with P>(r, r') = Y.lu x!o LlJLlL (r, r'; 0)m^ z (r'; 0). 
In matrix notation, Eq.[15]can be expressed as: 

rO = m z (16) 

which provides a means of calculating of U: 

0=r- 1 m z (17) 

Eq. [17] allows us to generate U through knowledge of only the ground state magnetization 
and the Kohn-Sham susceptibility Xo- An analysis of Eq. Q] shows that in the absence of 



an external magnetic field parallel to the z-direction the full dynamic susceptibility x will 
have a pole at zero frequency, if in fact U is generated from Eq. [17] Thus, by this scheme 
we generate an effective U compatible with the Goldstone theorem. Stated otherwise, the 
correct U is the one with the lowest eigenvalue of the denominator of EqfTJ associated with 
the magnetic moments as components of the eigenvectors. In the following we shall show 
through explicit calculation that the prescription in Eq. [17] can be applied to clusters of 
moment bearing ions which consists of dissimilar atoms. 

III. THE MASTER DYSON EQUATION WITHIN TD-DFT 

Let us briefly derive the master Dyson equation which leads to Eq. [TJ within the TD- 
DFT. By applying a linear variational approach, one assumes similar initial conditions as the 
ones in the previous section: i.e. a magnetic system with an initial charge density n (f), 
a magnetization pointing along the z-direction and an exciting time-dependent transverse 
magnetic field B ext (f; t) with small magnitude that allows us to use linear response theory. 
The result is an induced transverse magnetization localized in the (xy) plane perpendicular 
to the direction z. The art of TD-DFT is to relate and connect the induced transverse 
magnetization m xy (F; t) to the externally applied magnetic field. The dynamic susceptibility 
we seek may be expressed as a functional derivative of the transverse moment with respect 
to the external field, evaluated at zero external field: 

5m'[B ext ](rt) 



x IJ (rt,rt') 



(18) 

B e xt=0,rio 



5Si xt (Pt') 

where x is the response function we seek. In regard to the superscripts /, j and the definition 
of the vectors r, f see the remarks after Eq. El The convention we use here is the same as 
that employed for the single particle Green function. 

Within the atomic sphere approximation (ASA) and assuming once more an applied mag- 
netic field with spherical symmetry within the unit cell we may write 

<y(tt) = E/ dr j dt' X iJ (rt,rt')Bi xt (rt'), (19) 
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Upon resorting to the spherical harmonic expansion discussed above, this becomes 



m'Jrt) 



E 

j 



df / dt' 



LU;L 2 L 3 



X' LLl ' jL2i3 (rt, r't') x 



Y L {?)Y Ll {?)Y L2 {?)Y Lz {f)B ] ext {r't') 



(20) 



where r and r' are the magnitude of the vectors 7 and 7. 

If we integrate both sides of the previous equation over df we find: 



47r< y (rt) = J2 fdr' [ dt'J2x iLLlJUL (rt, r't')B{ xt {r't') (21) 

J J J LU 



Thus the functional derivative given by Eq.[TH] could be simplified to 

5m' xy [B ext ]{rt) 



X u (rt, r't') = 4tt- 



(22) 



B e xt=0,rio 



5B J ext (r't') 

where we define x u = Sll, x lLLljLlL - The same procedure is repeated for the magnetic 
response function Yn of the Kohn-Sham non interacting system which involves not only B ext 
but B e ff as well(14l]: As mentioned previously, B eff is the magnetic part of the effective 
Kohn-Sham potential (V eff — Vj ff ). After a Fourier transform with respect to time we obtain 
a form that maps our calculation onto the same structure employed many years ago by Lowde 
and Windsor 124J1. This remains often used in recent tight-binding simulations of magnetic 
excitations [6|, [7(] where it is found that the scheme accurately reproduces results found 
through use of a more sophisticated description of the Coulomb integrals. Our derivation 
elucidates how the structure introduced by Lowde and Windsor emerges from TD-DFT. 
We now have 

X! j {r,r';u) = x!i(r,r';u) 

+ X) / dr " [ dr '"Xo( r ' r ": u)U kl (r", r"'\ u)x! J {r"', r; u) (23) 

where the integrations are only over the magnitude of 7 and 7, with the site labeled ma- 
trix function shown. The effective Coulomb interaction U' J (r,r';u) may be expressed as a 
functional derivative given by 

ffS; ff (r;w) 



U IJ (r, r';uj) 



47T(5/TV (r / ; u) 



(24) 



B e xt=0,riQ 



Within ALDA prescription of the transverse response of the spin system, Eg. [241 simplifies 
to |27fl 



U ' J ( r > r '>") = f eff if''°n^ d IJ . (25) 
4tt m' z (r, 0) 

The object in Eq. 1251 will be noted as Udft is in the litterature often referred to as the 
exchange and correlation Kernel K xc . This is, it should be noted, exactly the form derived in 
Eq. [13] extracted from the sum rule Eq. [101 

From Eq. [25] it is obvious that U could be considered as a local exchange splitting divided 
by the magnetization. 

IV. CALCULATION OF THE KOHN-SHAM SUSCEPTIBILITY 

As shown in Eq. [21 the Kohn-Sham dynamical susceptibility is a convolution of two Green 
functions. The function Xo can be separated into a sum of two terms: l x which involves 
Green functions that are analytical in the same half complex plane, so / x itself is analytic, 
and then one has l 2 which is non analytic |6[|. For positive frequencies: 



F 



/?(f, r- u) = ±- I dzf(z) ( Gjj(K ?; z + u )Gj.(r, r; z) 



(26) 



- G}*(r>,r,z)Gj*(r,r>;z-u)\ 
and 

q{r,r;u) = y " dzf(z)(^-Gf J (r,r f ;z + Lj)G}*(r,7*;z) 

+ Gf i (r,f; Z )Gj*(r,r';z-u)] (27) 



Such a separation is attractive since / x can be calculated through use of a regular energy 
contour in the complex plane [28] with a modest k- and energy-mesh. In Ref. j^, Q], the 
energy contour consists of a line perpendicular to the real-axis starting at the Fermi energy 
and going to infinity. This is unfortunately not possible with the KKR-method since unwanted 
core states would then be included. Thus, the lower limit of the energy integration is chosen 
well below the valence band minimum. I 2 can be calculated along a line parallel to the real 
axis. This requires usually a very substantial numerical effort since a large number of k-points 
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as well as a dense energy mesh are needed. However, it can be shown that the integration 
is limited to a small energy controlled by u. In our discussion of spin excitations we are 
interested in frequencies u small compared to bandwidths, so the integrations involved in 
l 2 can be carried out readily. The computational effort is thus enormously reduced. Upon 
introducing a variable change we may write: 



F 



/«(?, f- u) = -— / dzGfjir, r;z + oj)GJ*{r, 7; z) (28) 
2^ Je f -w 

The use of two different contours can lead to a slightly different treatment of rather 
similar terms in l 1 and / 2 . In order to improve numerical stability, in the present analysis the 
two terms are arranged so they differ a bit from those presented in Ref. [gj. We write 



!'{{r,r-u) = — dz 



f{z)Gf J (r,r f ;z + uj)Gj l (r',r,z) 



f(z + u)Gf*(7, r;z + u)Gj*(K r<; z) 



ji 



(29) 



+^ J e u dzf{z)Gi(r,f , -,z + u)G] l {f,r,z) 



The second term on the right hand side of the previous equation can be added to l 2 which 
leads to 

7j'(r, f; u) = ^ jj ^ dzGfjir, f- z + u){Gl{7, r; z) - Gj*(r, f; z)) (30) 

while 

TiiKf-.u) = — J dz^f(z)Gf J (r,r-z + u)Gj l (r f ,r-z) 

-f(z + oj)Gf*(f*, r;z + oj)Gj*(r, 7; z) ) (31) 



or 



Ti(r,7;u) = ± p dz(f(z-u)Gi(r,f;z)Gj i (r,r,z-u) 



-f(z)G};(f, r; z)G];(T, f; z - u)j (32) 

This procedure just outlined is found to be stable and requires to calculate one less Green 
function. Up to now we have considered positive frequencies u. 
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Negative frequencies lead to slightly different forms of /i and / 2 : 

tiT.r-.u) = ^ p d Z f(z)Gf J (r,r;z-u)G] l (r,r,z) 

- f(z - u)Gf*(f, r; z - u)Gj*(K f; z)) (33) 

and 



T 2 {f, f; u) = ^-j dzGfjir, f; z - uXtftf, r; z) - Gl*(r, f; z)) (34) 



These expressions can be evaluated with modest numerical efforts since the required Green 
functions are the same than those calculated for the susceptibilities at positive frequencies. 

V. AN APPROXIMATE FORM FOR THE SINGLE PARTICLE GREEN FUNCTIONS 



The Green functions are provided by the KKR-GF method 26] 



G u {f,r-z) = J2- i ^R iL (r < : Z )H iL (r^z)5 u , LLl + R iL (n Z )G i B LJU (z)R^^: Z ) (35) 

where G B is the structural Green function. Here the regular R and irregular H solutions of 
the Schrodinger equation are energy dependent, and this makes the calculation of Xo in Eq. [TJ 
tedious and lengthy. Thus, instead of using Eq. [35] while evaluating Xo. we introduce the 
following simplification that captures the physics central to the systems of interest to us. 
In its spectral representation, the Green function is given by 

r ( ^ v CLjiE^Rjir, E^jE^Rj^F; Eg) 

k LU k 

where R' L {r, Eg) is a suitably normalized solution of the Schrodinger equation within the unit 
cell /'. 

Various Ansatz can be proposed to simplify the previous form. Instead of working with 
the energy dependent wave functions, one could use an energy linearized form of the wave 
function as done, for example, in the Linear Muffin Tin Orbital method [29(] or in the Full 



Potential Linearized Augmented Plane Waves method [30(]. Our Ansatz expresses the Green 
functions in terms of energy independent wave functions such that: 

G u {r, r;z) ~ — * (37) 

k LU k 
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or 

LL 



G U {r,f;z) ~ ^tiWG^Wfcin (38) 

with 



Note that after modifying the wave functions we naturally replaced the amplitude a by a 
different one (/3). 

Since our KKR-GF method generates the full Green function as given in Eq. [35l one could 
calculate G^ Ll (z) from 

IJ { 3 J dr(p' L *(r)(p' L {r) J dr'(pj L i(r>)4>i L i*(ri) { ' 

where on the right hand side of Eq. |40]we insert the full KKR Green function displayed in 
Eq.[35j 

The terms in the denominator are normalization factors. Thus, instead of working with 
<fi lL (r) we introduce 

*'H0 = ^ r (41) 



/ dr(p' L *{r)(p' L {r) 

where we choose 4>' L {r) = R' d (r;E F ), i.e., the c/-regular solution of the Schrodinger equation. 
This is appropriate for the calculation of the c/-block of the susceptibility. We propose here an 
expansion in terms of energy independent d like wave functions we choose to be the regular 
solutions of KKR-GF theory evaluated at the Fermi energy. Our focus is on low energy 
excitations of 3d moments so as we shall see below this choice is appropriate. 
Within the KKR-representation of the Green function G,y(z) is evaluated from: 



O z ) = E r d?H iL >{r>;zW L {r) f dr^ L ^(r)R' L ^r-z)5 IJ -L 2 L 3 (42) 

L 2 L 3 \ J ° J ° 

/rws pr W s 
d?R'L(7;zW L {r) J drilj iLl *(r)H iL z(nz)5 iU2L3 



drip iL *{r)R lL2 (r,z)G l B L2jL3 (z) j d?R JL3 (f;z)^ Ll (f)) 



where r ws stands for Wigner-Seitz radius. 
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VI. THE FINAL DYSON EQUATION 



Assuming the expansion in terms of energy independent wave functions described previ- 
ously, the final Dyson equation simplifies after some straightforward algebra into a strictly 
site dependent equation 

^ = + (43) 
where the d-block of the dynamical susceptibility is given by 

5g(r, r'; u) = ^V)<>)X>)<V)<V) (44) 

and 



piws 

TJ' = / drilt'firWftfU 1 {r)ip'i{r)ip\ d *{r) (45) 
Jo 

Within ALDA, we use Eq. [55] in Eq. [45] and obtain 

/r W s Di ( r n\ 

t/nftf (0 (46) 

If we want to use the sumrule we expand the susceptibility given in Eq. [T2] in terms of d-bloch 
susceptiblity expressed in Eq.|44]and repeat the same procedure used in section [Q] to find 

V = T~ 1 M Z (47) 

as written in matrix notation and r' J = x'o(0)M J z with M' z , calculated from the projection 
scheme proposed in section [V] is the magnetic moment of atom /'. U can be calculated 
once for every atom either from the previous sum rule, Eq. [47] or from Eq. [46] It can be 
understood as a Stoner parameter and gives once more a justification for the approach used 



by Lowde and Windsor 
by only one parameter. 



24]]: i.e. the effective intra-atomic Coulomb interaction is expressed 



VII. APPLICATION OF THE FORMALISM TO EXPLICIT EXAMPLES 
A. Single Adatoms 

We choose as an application of the formalism developed above the investigation of 3d 
adatoms and dimers positioned on the fourfold hollow sites of Cu(OOl) surface. In this 
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section, we focus on single adatoms. The calculations consist of the self-consistent de- 
termination of the electronic structure of these nanostructures using the usual KKR-GF 
scheme [26]. Once this is done, we generate the Green functions needed to calculate Xo, 
for the elements that bear a magnetic moments (Cr, Mn, Fe and Co), following Eq. [21 U is 
calculated either from Eq. [47] or Eq. [46] It is convenient to note that for the case of a single 
adatom /, Eq. [47] simplifies to U, = 4t at OJ = 0. 

We have already examined the spin dynamics of these systems in Ref. [211] where we have 
shown that the Green functions extracted from our approach (Eq. |43J) nicely reproduces the 
magnetic moment of the adatoms as calculated from a full DFT calculation. That this is so 
is illustrated in Fig. QJa). Indeed, interestingly, the ^-contribution to the total moment is, 
as expected, the most important and seems to be nicely reproduced by the projection of the 
Green functions into our choice of wave functions. 

We did not, however, discuss in Ref. [21J] the differences between values of U calculated 
from both schemes mentioned previously. In Fig. QJb) we show the values of U for the 
adatoms we have investigated. We find values of U very close to leV/ ji B for all cases we 



have studies. Himpsel 



251 ]. in his analysis of a large body of photoemission data on moment 



bearing 3d ions, has concluded that leV/ji B is a universal value that applies to diverse 
moment bearing 3d transition metal ions. As discussed in Ref. [211] . leV/ji B is also used 
commonly ETB calculations Q, Q]. Thus, we are pleased to see these values emerge from 
the scheme set forth here. The relative error or U values generated from density functional 
theory, as measured by the ratio ( UpFT ~ L,summle ) are depicted in Fig. HJc). The error is the 



DFT 



highest for Cr-adatom while the lowest is seen for Co. It is interesting that the observed 
error does not exceed 15% which is still much lower then what has been estimated by Sasioglu 
et a/.j^] while investigating bulk Ni. 

In Fig. [21 we show examples of the imaginary part of x for a Mn adatom positioned on the 
fourfold hollowsite of Cu(001) surface after applying an additional spatially uniform static 
magnetic field. The imaginary part of x describes on the resonant response of the local 
magnetic moment of Mn-adatom. As required by the Goldstone theorem, a zero frequency 
resonance is expected when no DC field is applied. We have verified numerically that this 
feature is present, when our method of determining U is employed. As soon as a DC field 
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FIG. 1: (a) Comparison between magnetic moments (in /ig) of adatoms calculated by the full 



with those calculated from the proposed projection scheme discussed in the text and 



2l| . Values of —U (eV//ie) are shown in (b) calculated from either from Eq. [57] or from 



KKR-GI 
in Ref. 

Eq. [46] while in the insert (c) we plot the percentage error defined as the difference between Uoft 
and L/ sumru | e divided by Udft- 
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pointing along the initial direction of the moment is applied, as discussed many years ago 
the local response of the moment displays a g shifted Zeeman resonance, broadened very 
substantially by decay of the coherent spin precession to particle hole pairs, whereas the 
total moment of the system precesses with g=2 and zero linewidth. Thus, experiments 
such as STM that are highly localized probes of the dynamic response of the moment see a 
qualitatively different response than very long wavelength probes such as microwave resonance 
or Brillouin light scattering. In the latter methods, both g shifts and linewidths have their 
origin only in terms in the system Hamiltonian that break spin rotation invariance. Examples 
are spin orbit effects, along with coupling of spins to lattice degrees of freedom. 

We see in Fig. [2] that the resonant frequency scales linearly with the applied DC field, 
as does the width of the structure in the local response of the moment. The width of the 
resonances is controlled by the local density of states [31], and is thus strongly influenced by 
the position of the d levels relative to the Fermi energy. 



B. Dimers of Identical Adatoms 

Let us turn to the case of dimers. We consider two identical adatoms each adsorbed in 
nearest neighbor four fold hollow sites on Cu(100). At such distances, their interaction is 
modest compared to energies which characterize the one electron properties of the system. 

In Fig. [31 we show effective values of U generated by different means of selecting this 
parameter. The one calculated with use of Eq. [461 refereed to as U D f T , is systematically 
smaller than that which follows from the sum rule in Eq. [471 We saw the same trend in our 
earlier discussion of single adatoms. Of course, if one employs U D ft in the calculation of the 
dynamic susceptibility the Goldstone theorem is not obeyed. We turn next to a discussion 
the two choices U + and U- that appear in Fig. [3] 

We discuss local dynamic susceptibilities x 11 . X 22 . X 12 ar) d X 21 - The superscripts refer 
to atomic sites where the atoms in the dimer are located. The response function x' J gives 
the response of the moment at site / in response to a spatially localized field applied to site 
j. So far, everywhere, upper cases were used for / and j site labels in the susceptibility. For 
the case considered in this section, where each atom in the dimer is identical and there is 
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— MnadatomonCu(OOl) 




Energy (meV) 

FIG. 2: Imaginary part of the transverse dynamical magnetic susceptibility for a Mn 
adatom/Cu(001) surface. After applying different DC magnetic fields, resonances are obtained 
and are shifted to higher frequencies by increasing the magnitude of the field. The corresponding 
Zeeman frequency with g=2 for the fields chosen are represented by the black circles. Thus the g 
shift is negative for this example. 

reflection symmetry through the midpoint of the line that connects their centers, we have 
X 11 = X 22 ar| d also x 21 = X 12 ; l n the next section we consider a dimer formed from two 
dissimilar atoms, so the equalities just stated do not hold. 

The Goldstone theorem requires that in the absence of an externally applied field (and 
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Cr Mn Fe Co 

Adatoms on Cu(OOl) surface 

FIG. 3: Different values of U obtained with different schemes for Cr, Mn, Fe and Co dimers 
deposited on Cu(OOl) surface. See the discussion in the text for the discussion of the various 
criteria for choosing U. 

in the absence of spin-orbit coupling) each element x' J must have a pole at zero frequency. 
This is insured if U is such that the determinant D formed from the matrix 1 — Uxo vanishes 
at zero frequency. For our dimer that consists of two identical atoms we have D = (1 — 
UXo 1 ) 2 + (^Xo 2 ) 2 - Upon setting D = 0, we encounter a difficulty. The criterion yields two 
acceptable values of U, U + = (Xc^ + Xo 2 ) -1 ar| d U- = (Xo 1 ~~ Xo 2 ) -1 - ' n Fig. El the red curve 
provides values of U+, for the ions we consider, and the blue curve U-. The two values of U 
determined by this criterion are quite close to each other, because on the electron volt scale 
the interaction energy between the two moments in the dimer is quite small, as noted above. 

One then must address which of the two choices for U discussed in the previous paragraph 
is the proper physical choice. To see this, we must refine our criterion. For the dimer with 



19 



two identical atoms, we can make a decision which value of U is the proper choice. If 
we consider the mode structure of the dimer, there is an acoustical mode wherein the two 
moments precess in phase, and an out of phase optical mode we shall discuss below. The 
Goldstone theorem requires the acoustical mode to have zero frequency. Thus, it is the 
function x a = X 11 + X 22 + X 12 + X 21 that also must have a pole at zero frequency, since 
this describes the response of the total moment of the dimer to a spatially uniform applied 
transverse field. For our simple dimer formed from two identical atoms, it is a simple exercise 
to find an expression for x a ■ O ne nas x a = (Xo 1 + Xo 2 )/[l ~~ ^(Xo 1 + Xo 2 )]- Thus for a pole 
to occur at zero frequency in this response function, we must choose U = U + . The sum rule 
provides us with the same criterion. 

For the case of the dimer just considered, it is straightforward to deduce the appropriate 
choice of U through examination of x a ■ However, for more complex arrays of spins the 
task of choosing U is not simple. Suppose, for instance we have N spins in the form of a 
one dimensional structure or possibly an island. From the numerical point of view, one may 
work with the analog of the determinant D discussed above. Exploration of its zeros at zero 
frequency will yield N possible values of U. Also if the spin structure consists of dissimilar 
atoms, each atom will be characterized by an appropriate value of U. As we shall see in 
the next section, the sum rule allows one to generate appropriate values of the interaction 
strength for each individual atom in a more complex structure. 

We turn next to the description of the spin dynamics of the dimer. For the dimer, 
we expect two resonances, an acoustical mode located obviously at u = and an optical 
mode at positive or negative frequencies. In general, the appearance of negative frequency 
modes in the dynamic susceptibility signal an instability of an assumed ground state. In 
the studies presented here, we assume a ferromagnetic ground state for the dimer. The 
appearance of a negative frequency optical mode is a signal that the atoms in the dimer 
are coupled antiferromagntically, so the ferromagnetic ground state is unstable. Thus the 
dynamic susceptibility can be used as a probe of local stability of assumed structures. 

It will be useful and interesting to compare our full dynamical calculations of the response 
of the dimer with the often used localized spin model, where effective exchange interactions 
are calculated within an adiabatic scheme. Such adiabatic scheme has already been used 
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for the investigation of different kind of systems (see e.g.Refs. |32h35|]). Through adiabatic 
rotations of the moments, [36], we extract an effective exchange magnetic interaction, J, by 
fitting the energy change to the Heisenberg form 

H = -Jei • e 2 (48) 

where and e 2 are unit vectors. By this criterion, we find that the ground state is antifer- 
romagnetic for Cr- (J = —19.8 meV) and Co-dimers (J = —14.9 meV) and ferromagnetic 
for Mn (J = 16.3 meV) and Fe (J = 30.4 meV). Since the dynamical susceptibility was eval- 
uated through use of ferromagnetic state for all the dimers, we expect an optical mode at 
positive frequencies for Mn and Fe dimers and at negative frequencies for Cr and Co dimers. 

We find that the dynamic susceptibility of the dimer is remarkably sensitive to the choice 
of the effective U. We see in Fig. [3]that numerically the difference between U + (= U sumru \ e ) 
and U- is quite small. Yet as illustrated in Fig. IUa), we show lm(x n ) calculated with the 
choice U = U-. For all four magnetic ions, the signature of the Goldstone mode is evident. 
For the Cr dimer, we see the clear signature of the optical mode at positive frequency. This 
suggests that, in contrast to the conclusion based on the adiabatic exchange analysis, the 
ferromagnetic ground state of Cr is stable. The optical modes of Mn, Fe all reside at negative 
frequency so for these three the results in Fig. 0Ja) suggest the ferromagnetic ground state 
is unstable. These results are also incompatible with the conclusions based on the adiabatic 
exchange integrals. 

In Fig. |4](b), we show results for lm(x n ) which follow from the choice U = U + . We now 
have results fully compatible with the conclusion based on the adiabatic exchange analysis. 
The sum rule has led to the correct selection of the effective U. 

Within the framework of the Heisenberg model, the optical mode should be an eigenmode 
of the system, and thus it will have zero linewidth. We see in Fig. HJb) that the optical mode 
for the Fe dimer and the Mn dimer have very substantial width. The origin of this broadening 
is in decay of the optical mode to Stone excitations. The itinerant character of the local 
moments is responsible for this linewidth, which elementary considerations suggest should 
increase linearly with the frequency of the optical mode. Thus, the linewidth of the optical 
mode of the Fe dimer is substantially broader than that of the Mn dimer. In the ground 
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state, hybridization between 3d states of the adatom and the conduction degrees of freedom 
on the Cu substrate results in "virtual levels" whose width is in the range of a few hundred 
meV. At the level of the spin dynamics, we see the large broadening of the optical mode as 
another reflection of the itinerant character of these systems. We note that in Spin-Polarized 
Electron Energy Loss Spectroscopy (SPEELS) studies of spin waves in ultrathin films very 
large linewidths are observed for high frequency, large wave vector modes [37]. The data is in 
excellent accord with theoretical calculations that assign the large linewidth to the damping 
by decay to Stoner excitations [7], very much as we see in the optical modes displayed in 
Fig. Kb). 

It is of interest to compare the frequency of the optical modes with the prediction of 
the Heisenberg model, with interspin exchange generated adiabatically as discussed above. 
If one considers two spin exchange coupled spins described by the Hamiltonian — J s Si • S 2 
the frequency of the optical mode is easily seen to be J s (Si + S 2 ). In Eq. |48] ei i2 are 
unit vectors, so J s = J/SiS 2 - Thus, in terms of the effective exchange couplings quoted 
above, with Si = S 2 = S the optical mode frequency is 2J/S. For the Mn and Fe dimers 
whose optical modes are illustrated in Fig.|4](b), the predicted frequencies are 15.4 meV and 
39.2 meV, respectively. The agreement with the optical mode of the Mn dimer is excellent, 
whereas the full dynamical calculation provides a somewhat smaller optical mode frequency 
for the Fe dimer. As discussed earlier, the coupling between the spin precession of the local 
moments and the Stoner excitations produces a mode softening not incorporated into the 
localized spin picture^, 7|. This coupling is considerably larger for the Fe dimer than the Mn 
dimer, as seen by a comparison of their linewidths. 



C. Dimers Formed from Different Adatoms 

We now turn our attention to a lower symmetry spin structure, dimers made of different 
magnetic adatoms. We study the MnFe dimer and the FeCo dimer, once again with the 
magnetic ions sitting in nearest neighbor fourfold hollow sites on the Cu(lll) surface. Here 
the two atoms do not have the same magnetic moments. Also the effective U is different for 
each atom. In this circumstance it is difficult to envision adjusting the values of U by hand 
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Frequency (meV) 

FIG. 4: local Imx 11 is shown for the four dimers based on: Cr-, Mn-, Fe-, Co- adatoms. To calculate 
X two possible schemes of evaluating are considered: either in (a) using U— or in (b) using 17+ . It 
turns out that corresponds to the value obtained from the sumrule (Eq. |4"T|) derived in the text 
while U that calculated from a simple iterative scheme out of Udft would converge to the wrong 
U when investigating Cr and Co dimers. The reason is that, for the latter elements, contrary to 
U+, U- is closer to Uoft- The optical modes, estimated for Mn and Fe from a Heisenberg model, 
are represented as dashed lines. 
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to obtain the zero frequency pole in the dynamic susceptibility. We have here a circumstance 
where the sum rule allows us to address the problem directly. Notice from Eq. |47]that though 
its use, we can determine the appropriate value of U for each atom in the dimer. Before we 
discuss imaginary part of the dynamical susceptibility let us discuss values of the magnetic 
moments and U's. 

TABLE I: Comparison between magnetic moments (in jLtg) and values of U's (eV//^e) for dimers 
made of different adatoms: MnFe- and FeCo dimers. 

Mn/Fe Fe/Co 
M d : projection model 3.85/2.74 2.78/1.64 

M to t a i 4.23/3.06 3.13/1.82 

-U DFT 0.89/0.94 0.94/0.95 

"^sumrule 0.97/0.98 0.98/0.98 



In Table [D the magnetic moments calculated with our projection scheme are shown and 
compared to the values that follow from the full KKR treatment of the ground state. In 
the first line of Tabled] the moment which appears is the contribution with d-\\ke symmetry, 
since this is the portion built into our Ansatz for the Green function used to compute the 
Khon-Sham susceptibility. It is interesting to note the substantial difference between the 
magnetic moments of two adatoms in the dimer. It is the case here as for the single adatom, 
the U calculated from Eq. |46] understimates the value of U needed to realize the Goldstone 
mode. From Eq. [471 we may deduce the value U, for each of the adatoms in the dimer. We 
find 

i o 90 

U, = m 1 ° — (49) 

■L 1 ? ?1 11?? v^-v 

Ao AO Ao Ao 

and 

m |Xo Xo 

2 19 oi 11 99 V / 

Ao Ao Ao Ao 

It is interesting that the sum rule gives similar values of U for both atoms in the dimer, and 
also that U is very close to leV/ ji B - That this is so is very compatible with the conclusion 
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of Ref. [25], which is based on an empirical study of photoemission data on 3d transition 
metal ions in diverse environments. 

The mapping to the previously defined Heisenberg model predicts a ferromagnetic ground 
state for both dimers investigated. Indeed the magnetic exchange interaction is positive in 
both cases with JMnFe = 28.1 meV (Heisenberg frequency 31.6 meV) and ^FeCo = 12.5 meV 
(Heisenberg frequency 21.7 meV). This indicates, as discussed above, that the imaginary 
part of the dynamical magnetic susceptibility for every adatom should show a resonance at 
positive frequencies that is the signature of the optical mode. In Fig. [5](a) and (b) we plot 
X 11 and x 22 respectively for the FeCo- and MnFe-dimer. 

A most striking feature of the results displayed in Fig. [5] is that the peak positions in x 11 
and x 22 occur at distinctly different frequencies. This is particularly clear in Fig. E£b), where 
the influence of damping is somewhat more modest than in Fig. [5fa). We see that the peak 
in x FeFe °ccurs at 30 meV, whereas that in ^ MnMn is distinctly downshifted to 27 meV. 

This behavior is at variance with the Heisenberg description of the excitation spectrum 
of two well defined localized spins. As we have seen, if we have two well defined, localized 
spins coupled together by the exchange interaction — J<?i ■ e 2 , the pair has two excited states 
associated with small amplitude motions, the acoustical mode at zero frequency (which we 
see in Fig. [5]) and the optical mode at the frequency 2J/S. Thus, the optical mode peak 
in the excitation spectrum for each member of the dimer should be at exactly the same 
frequency, in this picture. While the oscillator strength of each peak will differ, there is a 
unique excited state energy of the pair. 

The shift in the peak positions evident in Fig. [5] is a consequence of the itinerant nature of 
the magnetic moments. As each moment precesses, as we have seen, the motion is damped 
heavily by the coupling of the moment to the Stoner excitations of the paramagnetic host. 
In the case of the FeMn dimer, the motions of the Fe spin are damped far more heavily 

n 

that those of the Mn spin, as we may appreciated from Fig. 1(b) of Ref. |21[|. This has the 
consequence that the peak in lmx MnMn is dragged down to a frequency somewhat lower than 
that in lmx FeFe - We may see this by constructing a toy model that consists of two Heisenberg 
coupled spins, each of which is coupled to a reservoir that produces damping a of the form 
encountered in the Landau-Lifschitz-Gilbert equation. The linearized equations of motion for 
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this system reproduces the offset in the peaks evident in Fig.[5fb). We illustrate this in Fig. [6] 
where Imx 11 and Imx 22 mimic the imaginary parts of ^ MnMn and x FeFe - By increasing the 
strength of the damping parameter a 2 compared to a±, we observe a shift to lower energies 
of the optical mode in Imx 22 {i.e. lmx MnMn )- it is striking to observe the completely different 
shape of the optical mode of Mn-spin just by modifying a neighbor. Indeed, by comparing 
the optical mode observed in | m ^- MnMn we observe also that it is much more heavily damped 
in the mixed dimer MnFe (Fig. E^b)) than in the pure MnMn dimer (Fig.@Jb)). The physical 
reason behind this intriguing behavior is that in the MnFe configuration, the Mn-spin during 
its precession feels the magnetic force of the heavily damped Fe-spin which provides more 
damping on Mn. It would be of great interest to employ STM based spectroscopy to explore 
the response of the two spins in a dissimilar dimer such as that just discussed. 

VIII. CONCLUSION 

We have developed and presented a theory based on TD-DFT and the KKR-GF method 
to extract dynamics magnetic susceptibilities of moment bearing adatoms and adatom dimers 
on surfaces. In our method, the electronic structure is described within an ab-initio scheme 
with KKR Green functions as the basis. Thus, no parameters need to be introduced, as 
in studies that employ the empirical tight-binding method. As important feature of our 
approach is that it may be implemented with a modest expenditure of computational effort. 
It is thus suitable for exploration of complex magnetic structures on surfaces that contain 
several magnetic ions. In this paper, we illustrate the method with application to magnetic 
dimers formed from either identical or dissimilar adatoms. 

As discussed above, a difficulty with past TD-DFT studies of spin excitations not only 
on surfaces, but in bulk materials as well is that the effective value of the Hubbard U which 
emerges from the standard approaches is not compatible with the Goldstone theorem that 
guarantees that the low lying acoustical spin-excitation has zero frequency. This difficulty has 
led others to make ad-hoc adjustments in the value of U. A feature of the present analysis 
is the introduction of a sum rule from which proper values of this parameter emerge. This 
eliminates the need for ad-hoc adjustments. It should be remarked that in simple systems, 
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where the analysis can be phrased in terms of a single value of the effective U, it is not difficult 
to insure satisfaction of the Goldstone theorem through an ad-hoc correction, though in our 
view this is an unsatisfactory procedure that compromise the theory at the fundamental 
level. Additionally, for a multicomponent system, the ad-hoc correction procedure becomes 
problematic in practice. As we see from our discussion of the dimer constructed from two 
different magnetic ions, our sum rule approach is readily and easily implemented for multi- 
component systems. 
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Appendix 

In this appendix we provide a derivation of the useful identity presented in Eq. [6] 

The Green function G(z) of a Hamiltonian operator H is defined by the operator equation 

g = 7±h < 51 > 

If no spin-orbit coupling and non-collinear magnetism are considered, the previous equation 
holds for every spin-channel (f or \). Thus 

GtW = - (52) 

In addition we have: 

z - H l = z-rf + H i -H 1 ~ (53) 
that can be multiplyed from both sides from the left by (z — H^) -1 and from the right by 



(z - H 1 ") -1 . This leads to 



1 1 + ^TJT(^ - H^)-l- (54) 



z - Ht z - H± z - H± y J z - Ht 
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i.e. 

G t = G^ + G^BeffG^ 

where we define B eff = H l - /-/ f . 
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FIG. 5: local Imx for dimers with mixed adatoms are shown in (a) for FeCo dimer and in (b) 
for MnFe dimer. Eq. 07] based on the sum rule derived in the text was used to define U. It is 
interesting to note the presence of resonances at positive frequencies expressing a ferromagnetic 
ground state for both dimers. Within each dimer, the pics related to every adatom are not located 
at the same position since the g-shift depends on the nature of the adatom. 
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FIG. 6: The response function Im(x 1:L ) and Im(x 22 ) for two spins of unit length coupled by an 
exchange interaction of strength J = 1. Here, we mimic Fe and Mn by considering each spin 
coupled to a reservoir that provides a damping parameter oci^ (1 for Mn and 2 for Fe) whose 
values are given in the inset. 
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